Summary

This document guides users through the steps needed to identify calving grounds using a population-based method that we developed for barren-ground caribou (Rangifer tarandus groanlandicus) in the Qamanirjuaq Herd (Couriot et al. In prep). Here, we apply this method the Bathurst caribou movement data, available to registered collaborators on MoveBank (https://www.movebank.org). The method takes advantage of the observation that herd-level ranging areas and movement rates (speeds) tend to have two fairly distinct break points around the calving season: (I) during spring migration the rate of movement of the herd is high and the overall “ranging area” is high as there is consierable variation in the arrival times of individuals, (II) during calving both herd-level measures are quite small, and (III) post-calving, both speeds and ranging areas increase again, but neither typically to the extent of the migration season. In summary, the steps for implementing the method are to:

1. load and install the TuktuTools package;
2. download the movement data of Bathurst caribou on Movebank;
3. process the data by limiting the dates of analysis, removing outliers, and computing the daily mean locations; 4. estimate the daily movement rate of the herd (i.e. the average movement rate of all the females monitored at a given day) and the daily ranging area of the herd (i.e. the area used by all the females monitored during a given day);
4. identify the calving and early rearing period by segmenting the bivariate time series of the herd; and
5. identifying the calving grounds by estimating the spatial range of the herd during that calving period.

Background

Migratory tundra caribou, or barren-ground caribou, belong to sub-populations called herds which are defined by their traditional fidelity to specific calving grounds. These calving grounds are reached by individuals from widely dispersed wintering grounds during the spring migration season. Once on these calving grounds, there is a high degree of synchrony in calving, with 50% of the animals in a given herd giving birth within 4-7 days (Lent 1974; Couriot et al. In press). This behavior likely evolved as a strategy to mitigate predation, via remoteness, spatial aggregation and synchrony of births, while also maintaining access to forage-rich post-calving and summer grounds (Bergerud 1996; Cameron et al. 2020; Gunn and Miller 1986). There is evidence to suggest that caribou cows show memory-based fidelity to regions known to meet their energetic requirements for sufficient forage and predator evasion during these critical periods (Cameron et al. 2020). While caribou calves can stand on their legs 30 minutes after birth, they cannot keep pace with adults for several days (Lent 1974). The early post-calving period is also a critical period during which calf-cow pairs have high nutritional requirements (Lent 1974) and calves nurse frequently – up to 30 times a day (Russell, Kofinas, and Griffith 2002). This early post-calving period lasts approximately three weeks until the recovery phase, when calves start to eat green forage (Lent 1974; Russell, Kofinas, and Griffith 2002). For all of these reasons, movement rates of mothers are considerably lower during this early rearing stage. Once calves are fully mobile - typically within a few weeks of birth - females regain movement rates similar to those before calving. At this stage, animals are often highly aggregate as they coalesce into movement corridors towards insect relief or summer rearing areas (Lent 1966).

Given the social behavior and the low movement rate of female caribou during calving and the early-rearing period, we found it is possible to identify the calving and early calf-rearing period by simultaneously analyzing both the aggregation of females as well as their movement rate through time. Before calving, females are widely dispersed and moving at a high rate during spring migration, and the highly aggregated and slow-moving during calving and early rearing, and still fairly aggregated but moving much more rapidly after leaving the calving grounds as they depart the calving grounds for insect-relief and summer foraging areas. There can be considerable variability in these metrics, as animals may be more or less dispersed during migration and, in particular, post-calving, and in some years not animals make it to the calving ground to calve (Couriot et al. In press), but with a sufficient number of tracked animals, these periods can usually be distinctly separated along these two dimensions.

It is important to note that we define the “annual calving period” and “annual calving ground” to be the ecologically relevant interval of time and spatial area used both for calving and early rearing.

Our population-based method for determining the calving (Couriot et al. In prep) segmenting a joint (bivariate) time series of daily mean herd movement rate and ranging area. The daily movement rates of the herd at a given year is calculated by taking the individual displacement between subsequent daily centroids (i.e., the centroid of all the locations for a given female at a given day) and averaging those displacements across all females for each day between 15 May and 15 July. We expect the herd movement rates to show a pronounced dip as animals sequentially began to calve, followed by a slower ramping up of movement rates to levels near pre-calving. The daily ranging area is an estimate of the total spatial extent of the herd at a given day, which varies greatly throughout the year: it is greatest during spring migration, smaller as the animals aggregate on the calving ground, and smallest as the animals form post-calving aggregations as a mean to decrease insect harassment. To estimate the daily ranging area of the herd, computed the 90% contour of the Local Convex Hull (LoCoH, (Getz et al. 2007), using the tlocoh package in R, which has higher convergence properties than parametric kernel methods (Getz et al. 2007)). The daily LoCoH is the union of all the convex hulls fitted around each daily centroid of the females and 8 nearest collared female neighbors, on a given day. For example, for a given day in 2020 (where there were n=36 collared Bathurst females), the daily LoCoH was the union of the 36 convex hulls associated with each daily centroid and its 8 nearest neighbors. We chose 8 nearest neighbors as the best compromise between the lowest number of daily centroids per day (i.e., 13 females monitored in 2013) and the definition of the resulting polygon (Getz et al. 2007). The bivariate time series is then segmented using a two-dimensional time-space clustering method, implemented in the segclust2d R package (Patin et al. 2020). The estimated break points are then used to delineate the start- and end-dates of the calving season (Fig. 1).

Fig. 1 Figure illustrating the bivariate segmentation of movement rate (displacement per day) and square root of ranging area (km²) into three phases: pre-calving (spring migration, blue), calving and early post-calving (EPCP, red), and post calving aggregation (green), for the Qamanirjuaq Herd in 2018 (Couriot et al. In prep).

Once the herd-level estimates of the start and end time of the calving and early rearing period are determined, we define the greater annual calving areas and central annual calving areas as the 95% and 50% kernel utilization distribution (UD), respectively. These are estimated using the first and last location of the estimated calving season for each female for a given year. The set of locations used to determine these areas is includes the first and last location for each female during the season to minimize potential issues with over determined areas - an issue that can arise if all the movement data is used within the period due to high auto-correlation in the data.

Implementation

devtools::install_github("ocouriot/TuktuTools", build_vignettes = TRUE)

# Windows:
install.packages("tlocoh", dependencies=TRUE, 
repos=c("http://R-Forge.R-project.org", "http://cran.cnr.berkeley.edu"))

# Mac:
install.packages("tlocoh", dependencies=TRUE, repos=c("http://R-Forge.R-project.org", "http://cran.cnr.berkeley.edu"), type="source")

require(TuktuTools)

require(move)
mylogin <- movebankLogin(username='XXX', password='XXX')

bathurst.move <- getMovebankData(study = "GNWT North Slave Barren Ground Caribou: Bathurst",
    login = mylogin, removeDuplicatedTimestamps=TRUE)

myCRS <- "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=65 +lon_0=-120 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
bathurst <- process_moveNWT(bathurst.move, mylogin = mylogin, CRS = myCRS)

Simple feature collection with 1102642 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -204698.8 ymin: -519079.6 xmax: 1095812 ymax: 653533.9
CRS:           +proj=lcc +lat_1=50 +lat_2=70 +lat_0=65 +lon_0=-120 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
First 10 features:
     ID nickname                Time      Lon    Lat sex Year study_site                   geometry
1  X100    BA100 1996-04-11 11:46:00 -109.785 64.202   f 1996   Bathurst POINT (487700.3 -49908.54)
2  X100    BA100 1996-04-18 12:11:00 -109.744 64.429   f 1996   Bathurst  POINT (485765.5 -24931.5)
3  X100    BA100 1996-04-25 12:14:00 -109.525 64.255   f 1996   Bathurst POINT (499085.9 -42200.59)
4  X100    BA100 1996-05-02 11:20:00 -109.104 64.804   f 1996   Bathurst  POINT (509006.8 20650.65)
5  X100    BA100 1996-05-09 10:30:00 -108.937 66.153   f 1996   Bathurst  POINT (491856.1 168552.9)
6  X100    BA100 1996-05-16 11:19:00 -109.086 66.855   f 1996   Bathurst  POINT (472493.3 243917.4)
7  X100    BA100 1996-05-23 10:53:00 -109.201 66.835   f 1996   Bathurst  POINT (467920.9 240915.1)
8  X100    BA100 1996-05-29 13:12:00 -109.626 66.895   f 1996   Bathurst  POINT (448620.3 244502.2)
9  X100    BA100 1996-05-30 11:13:00 -109.373 66.809   f 1996   Bathurst  POINT (460995.4 236865.5)
10 X100    BA100 1996-05-31 10:49:00 -109.735 66.804   f 1996   Bathurst  POINT (445507.7 233826.7)
          x         y
1  487700.3 -49908.54
2  485765.5 -24931.50
3  499085.9 -42200.59
4  509006.8  20650.65
5  491856.1 168552.85
6  472493.3 243917.38
7  467920.9 240915.14
8  448620.3 244502.18
9  460995.4 236865.54
10 445507.7 233826.69

# Filter females and movement data between mid May and mid July from 2010 to the most recent data
bathurst_clipped <- bathurst %>% subset(study_site == "Bathurst" & sex == "f" &
                                  Year >= 2010 &
                                  ((mday(Time) >= 15 & month(Time) == 5) |
                                  (mday(Time) <= 15 & month(Time) == 7) |
                                  month(Time) == 6)) %>% mutate(yday = yday(Time))

# Filter females and movement data between mid May and mid July from 2010 to the most recent data
bathurst_clipped <- bathurst %>% 
  subset(study_site == "Bathurst" & sex == "f" & Year >= 2010) %>%
  prepData(start = "05-15", end = "07-16")

## Period clipped to 05-15 - 07-16
## Number of excluded individuals-years: 79

# example
scan_tracks(bathurst_clipped %>% subset(Year == 2018),  legend = TRUE, ncol = 3, cex = 0.6)

bathurst_filtered <- bathurst_clipped %>%
  subset(!(ID %in% c("BG236") & Year == 2012)) %>% 
  subset(!(ID %in% c("BGCA18129") & Year == 2018)) 

bathurst_dailymean <- bathurst_filtered %>% getDailyMean %>% 
  st_as_sf(coords=c("x", "y"), crs = myCRS) %>%
  mutate(x = st_coordinates(.)[,1], y = st_coordinates(.)[,2])

bathurst_herdmovements <- bathurst_dailymean %>% as.data.frame %>% 
  mutate(Date = as.Date(Time)) %>% 
  ddply(.(Year), getSpeed) %>% group_by(Year, yday, Date) %>% 
  summarize(speed = mean(speed, na.rm = TRUE)) %>%
  as.data.frame

head(bathurst_herdmovements)

bathurst_dailyarea <- bathurst_dailymean %>% ddply("Year",
  getDailyRange, crs = myCRS, method="LoCoH", nn = 8, level = .9) %>%
  st_as_sf(crs = myCRS)

head(bathurst_dailyarea)

daily_area_movement <- 
  merge(bathurst_dailyarea %>% as.data.frame %>% mutate(geometry = NULL),
        bathurst_herdmovements, by = c("Year", "yday")) %>% 
  mutate(area = area/1e6) %>% na.omit

calving_period <- daily_area_movement %>% 
  ddply("Year", getCalvingPeriod, min.duration = 5, 
        drawplot = TRUE, ylim = c(0,300))

calving_start_end <- bathurst_dailymean %>% 
  base::merge(calving_period[,c('Year','start','end')], by = "Year") %>% 
  subset(yday == start | yday == end) %>% arrange(ID, Year, yday)

require(adehabitatHR)
CRS = st_crs(bathurst_dailymean)

# 95% kernel density
greater_CalvingGrounds <- calving_start_end %>% group_by(Year) %>% 
  group_modify(~ getKernelUD(., h = "href", percent = 95)) %>% st_as_sf(crs = CRS) %>% 
  st_transform(crs=4326) %>% mutate(Year = as.factor(as.character(Year)))

# 50% kernel density
central_CalvingGrounds <- calving_start_end %>% group_by(Year) %>% 
  group_modify(~ getKernelUD(., h = "href", percent = 50)) %>% st_as_sf(crs = CRS) %>% 
  st_transform(crs=4326) %>% mutate(Year = as.factor(as.character(Year)))

require(mapview)
mapview(greater_CalvingGrounds, zcol = "Year")

st_write(greater_CalvingGrounds, dsn = "greater_CalvingGrounds.shp")
st_write(central_CalvingGrounds, dsn = "central_CalvingGrounds.shp")

# turn off the Spherical geometry (s2) to compute the union
sf_use_s2(FALSE) 

# get the union
greater_CalvingGround_pooled <- st_union(greater_CalvingGrounds)
central_CalvingGround_pooled <- st_union(central_CalvingGrounds)

# map the calving ground
require(mapview)
m1 <- mapview(greater_CalvingGround_pooled)
mapview(central_CalvingGround_pooled, map = m1@map)

Caveats and comments

Although the method we describe is this document has proven several benefits, it is worth noting some caveats. The method relies on the segmentation of the daily ranging area and the daily movement rate of the herd, which was easier and faster to compute than individual-based approaches and used readily-available computational tools in R. It aggregates the behavior of all individuals, and, in that way, points to some of the collective behavior of barren-ground caribou (Dalziel et al. 2015; Skoog 1968). However, although spring migrations are considered a social behavior, and occur relatively synchronously within a herd (Gurarie et al. 2019), females usually move in several groups with non-pregnant females migrating later and slower than calving females (Adamczewski, Gunn and Campbell, personal comments), which could lead to some uncertainty in the identified calving period’s dates. It is also worth noting that the daily ranging area of the herd is sensitive to the number of nearest neighbors chosen when computing the Local Convex Hulls (LoCoH). We recommend trying different number of nearest neighbors to compute the daily ranging area before doing the bivariate segmentation. Finally, this method has been developed to identify annual calving periods and calving grounds, to further be able to analyze inter-annual differences in size and location (Couriot et al. In prep). The multi-annual calving ground, consisting of the combination of all the annual calving grounds, is thereby sensitive to annual calving grounds location and area, which are known to vary inter-annually (Couriot et al. In prep). Given the different points stated above and since herds’ behaviors can vary inter-annually and among herds (migration timing, birth timing), due to variation in local environmental conditions Couriot et al. (In press), we recommend to visually inspect the resulting annual calving periods and calving grounds and to try several parameters to choose the ones that are the most appropriate to each herd. (e.g., the dates of the extended calving periods to filter the data before running the analysis, the number of nearest neighbors used to compute the Local Convex Hulls).

Acknowledgments

This method was originally developed in the context of a project entitled Analysis of barren-ground caribou calving locations in northern Canada, the goal of which was to define and identify the calving grounds of Qamanirjuaq herd caribou in Nunavut. This was a collaboration with Environment Climate Change Canada and supported by a grant from WWF. We are grateful to our collaborators on that project, Mathieu LeBlond (ECCC), Mitch Campbell (Government of Nunavut) and Olivia Crosby, all of whom were instrumental in developing this tool.

References

Supplementary materials

In order to simplify the reproduction of this tutorial to identify the calving period and annual calving ground for the Bathurst Herd, we provide hereafter an example of code for one year (here 2010) with minimal annotations.

# Load the packages
require(TuktuTools)
require(move)
require(adehabitatHR)
require(mapview)

# create MoveBank login  change username and password accordingly
mylogin <- movebankLogin(username='Ophelie Couriot', password='Rangifer')

# download data for the Bathurst Herd from MoveBank
bathurst.move <- getMovebankData(study = "GNWT North Slave Barren Ground Caribou: Bathurst",
    login = mylogin, removeDuplicatedTimestamps=TRUE)

# set the CRS
myCRS <- "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=65 +lon_0=-120 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

# process move data
bathurst <- process_moveNWT(bathurst.move, mylogin = mylogin, CRS = myCRS)

# Filter females and movement data between mid May and mid July in 2010
bathurst_clipped <- bathurst %>% 
  subset(study_site == "Bathurst" & sex == "f" & Year == 2011) %>%
  prepData(start = "05-15", end = "07-16")

# estimate the daily mean location for each individual
bathurst_dailymean <- bathurst_clipped %>% getDailyMean %>% 
  st_as_sf(coords=c("x", "y"), crs = myCRS) %>%
  mutate(x = st_coordinates(.)[,1], y = st_coordinates(.)[,2])

# estimate the herd's daily movement rate
bathurst_herdmovements <- bathurst_dailymean %>% as.data.frame %>% 
  mutate(Date = as.Date(Time)) %>% 
  ddply(.(Year), getSpeed) %>% group_by(Year, yday, Date) %>% 
  summarize(speed = mean(speed, na.rm = TRUE)) %>%
  as.data.frame

# estimate the herd's daily ranging area
bathurst_dailyarea <- bathurst_dailymean %>% getDailyRange(crs = myCRS, method="LoCoH", nn = 8, level = .9) %>%
  st_as_sf(crs = myCRS)

# create the bi-variate df
daily_area_movement <- 
  merge(bathurst_dailyarea %>% as.data.frame %>% mutate(geometry = NULL),
        bathurst_herdmovements, by = c("Year", "yday")) %>% 
  mutate(area = area/1e6) %>% na.omit

# do the bi-variate segmentation
calving_period <- daily_area_movement %>% 
  getCalvingPeriod(min.duration = 5, 
        drawplot = TRUE, ylim = c(0,300)) %>% 
  mutate(Year = year(Date))

# get the estimated calving period in 2010
calving_start_end <- bathurst_dailymean %>% 
  base::merge(calving_period[,c('Year','start','end')], by = "Year") %>% 
  subset(yday == start | yday == end) %>% arrange(ID, Year, yday)

# set CRS
CRS = st_crs(bathurst_dailymean)

# compute the calving ground in 2010 using the 95% kernel density
greater_CalvingGround <- calving_start_end %>% group_by(Year) %>% 
  group_modify(~ getKernelUD(., h = "href", percent = 95)) %>% st_as_sf(crs = CRS) %>% 
  st_transform(crs=4326) %>% mutate(Year = as.factor(as.character(Year)))

# compute the calving ground in 2010 using the 50% kernel density
central_CalvingGround <- calving_start_end %>% group_by(Year) %>% 
  group_modify(~ getKernelUD(., h = "href", percent = 50)) %>% st_as_sf(crs = CRS) %>% 
  st_transform(crs=4326) %>% mutate(Year = as.factor(as.character(Year)))

# visualize the greater calving ground in 2010
mapview(greater_CalvingGround)